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ABSTRACT 

Millisecond and binary pulsars are the most stable natural frequency standards which 
admits to introduce modified versions of universal and ephemeris time scales based 
correspondingly on the intrinsic rotation of pulsar and on its orbital motion around 
barycenter of a binary system. Measured stability of these time scales depends on 
numerous physical phenomena which affect rotational and orbital motion of the pulsar 
and observer on the Earth, perturb propagation of electromagnetic pulses from pulsar 
to the observer and bring about random fluctuations in the rate of atomic clock used 
as a primary time reference in timing observations. On the long time intervals the main 
reason for the instability of the pulsar time scales is the presence of correlated, low- 
frequency timing noise in residuals of times of arrivals (TOA) of pulses from the pulsar 
which has both astrophysical and geophysical origin. Hence, the timing noise can carry 
out the important physical information about interstellar medium, interior structure 
of the pulsar, stochastic gravitational waves coming from the early universe, etc. Each 
specific type of the low-frequency noise can be described in framework of power law 
spectrum model. Although the data processing of pulsar timing observations in time 
domain seems to be the most imformative it is significantly important to know to which 
spectral bands single and binary pulsars, considered as detectors of the low- frequency 
noise signal, are the most sensitive. Solution of this problem may be reached only if 
a parallel processing of timing data in frequency domain is fulfilled. This requires a 
development of the Fourier analysis technique for an adequate interpretation of data 
contaminated by the correlated noise with a singular spectrum. The given problem is 
examined in the present article. 

Key words: methods: data analysis - methods: statistical - pulsars: general, binary 



1 INTRODUCTION 

Millisecond and binary pulsars are known as exellent probes for testing theory of general relativity 
Taylor & Weisberg 1982, 1989), structure of interstellar medium (Rickett 1990, 1996) and interior 
of neutron stars (Cordes & Greenstein 1981, Kaspi et al. 1994) as well as setting upper limit on 
the energy density of primordial gravitational radiation (Kaspi , Thorsett & Dewey 1996, McHugh 
et al. 1997, Kopeikin 1997a, Kopeikin & Wex 1999). Rotational motion of a pulsar around its own 
axis had been proposed (Shabanova et al. 1979, Backer et al. 1982, Rawley et al. 1987, Matsakis 
et al. 1997) for using as a new time reference being analogue of universal time in astrometry. 
Quite recently, a new step toward to estiblishing a stable time scale on extremely long intervals 
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approaching 50-100 years has been suggested (Ilyasov et al. 1998, Kopeikin 1999). It is extracted 
from the orbital motion of pulsar in a binary system and represents the analogue of ephemeris time 
of classical astronomy introduced by Newcomb (1898)at the end of last century and based on the 
orbital motion of the Earth around the Sun. 

An adequate analysis of timing data requires deeper understanding of nature of a noise process 
dominating in pulsar timing residuals. As soon as the autocovariance function of the noise process 
is known the analysis in time domain becomes possible. Time domain analysis is the most informa- 
tive since the observed stochastic process is not usually stationary and includes the non-stationary 
component as well (Groth 1975, Cordes 1978, 1980; Kopeikin 1997b). Because pulsar timing obser- 
vations are conducted on relatively long time intervals the white noise of errors in measuring TOA 
of pulsar's pulses will be suppressed by the presence of a number of correlated, low-frequency (red) 
noises having different spectra and intencities. Henceforth, we are mainly interested in analysing 
the red noise. 

The simple model of such noise has been worked out by one of us (Kopeikin 1997b). It is based 
on the shot noise approximation and do includes a dependence of the autocovariance function on 
both stationary and non-stationary components of red noise. In the process of elaborating the given 
model a rather remarkable fact has been established (Kopeikin 1999), namely, that timing residuals 
and variances of some spin-down and all orbital parameters are not affected by the non-stationary 
component of the red noise at all. This discovery put on a firm ground the Fourier analysis of TOA 
residuals and variances of fitting parameters in frequency domain. This analysis gives exhaustive 
information about the noise process itself and admits us to reveal to which frequency harmonics in 
the spoectral expansion of the stochastic process pulsar timing observations are the most sensitive. 
Moreover, just we have an adequate approach for the treatment of red noise in frequency domain 
a lot of interesting applications is opened having the goal to study the physical nature of the low 
frequency noises with arbitrary spectrum. 

Any low-frequency noise can be approximately characterized by the power-law spectrum S(f) ~ 
f~ n where the spectral (integer) index n > 1. It is obvious that the spectrum has a singularity 
at zero frequency. Hence, the energy of TOA residuals comprised in such noise should has infinite 
value because the integral over all frequencies from zero to infinity taken from S(f) is divergent. 
Clearly, this has no physical meaning and one has to resort to special mathematical tricks in order 
to avoid this artificial divergency. There are two ways for curing this flaw. First of them consists in 
analytical continuation of the spectrum by changing it from S(f) ~ f~ n to S(f, A) ~ f~( n + A ) where 
A is the pure complex parameter being differnt from zero. Such model of the analitically continued 
spectrum gives convergent integrals which coincide everywhere on real axis with ones taken from 
the spectrum S(f) ~ f~ n exept for the point / = 0. In order to prescribe a physical meaning to 
such integrals we have to expand them in the Laurent series with respect to the parameter A and 
to take the finite part of the expansion. Such procedure has been used, in particular, by Kopeikin 
(1997a) for calculation of autocovariance function of stochastic noise of the primordial gravitational 
wave background having spectrum S(f) ~ f~ 5 (Mashhoon 1982, 1985, Mashhoon & Seitz 1991, 
Bertotti et al. 1983). 

The procedure of analytical continuation of divergent inetgrals is mathematically rigorous and 
theoretically powerful tool which gives well-defined and self-consistent results (Gel'fand & Shilov 
1964). However, it is inconvinuent for people doing numerical computations. For them, the second 
method of regularization of singular spectra seems to be more preferable and practically useful. 
It is based on the truncation of all divergent integrals at the lower cut-off frequency f — s along 
with corresponding modification of the power-law model of the red noise spectrum in order to 
avoid the model dependence of results of fitting procedure on the artificially introduced cut-off 
frequency. It will be shown that the simplest way of modification of the spectrum can be reduced 
to the addition to the existing power-law spectrum of red noise of the infinite sum consisting of the 
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Dirac delta function and its derivatives having local support at the lower cut-off frequency / = e. 
Such modernization of the spectrum preserves the structure of the autocovariance function and, 
as a consequence, do not vilolate results of numerical computations in time domain. This second 
method of regulariazation of divergent integrals will be used in the present paper. 

In what follows it is more convinuent for us to work in terms of dimensionless frequency and 
time. For instance, in binary pulsars it is preferable for analytic calculations to measure time in 
units of orbital frequency n& = 2n / P&, where P& is the orbital period of the binary. Then, frequency 
/ is measured in units of 1/P& and dimensionless time is the pulsar's mean anomaly u = n^r. 
Hereafter, we use u instead of r. 



2 REGULARIZED SPECTRUM OF LOW-FREQUENCY NOISE 

Any gaussian low-frequency noise is completely characterized by the autocovariance function which 
describes correlation between two values of stochastic process separated by the arbitrary time 
interval r = £2 — t\. Autocovariance function consists usually of two algebraically independent 
parts characterizing separately stationary R~(t) and non-stationary R + (t) components of the 
noise. Complete expressions of autocovariance functions for different events of low-frequency noise 
have been derived in the paper (Kopeikin 1997b) where the shot noise approximation of stochastic 
process has been used. Although both stationary and non- stationary components of autocovariance 
function are important for an adequate treatment of observations (Kopeikin 1999) we are dealing in 
the present paper only with the stationary part which posseses to be transformed into the spectral 
density of noise S(f) by means of the Wiener-Khintchine theorem 

rod 

R~(u)=2 S(f) cos(2tt fu)df. (1) 
Jo 

If (constant) intensity of noise is denoted by h n then autocovariance function of low-frequency 
correlated noise is determined by the expression (Kopeikin 1997b) 

{C n h n \u\ n ~ l , n = 2, 4, 6, random walk noise 
(2) 
C n h n u n 1 In |w| , n = 1, 3, 5, flicker noise 

where C n is a numerical constant of normalization. 

Functions, which might be appropriate candidates for the spectrum of noise procesess with the 
foregoing autocovariance functions, are S(f) = const/ - ™ where n is integer. However, integrals ([!]) 
from such power-law functions are divergent because of non-physical singularity at zero frequency. 
For this reason, regularization technique should be used because we don't know usually the low- 
frequency behavior of the spectrum. 



2.1 Analaitic Continuation Technique 

Analytic continuation regularization procedure is contained in that one extends the spectrum S(f) 
in the complex plane domain by introducing the function 

S(f,A) = constf- n - A (3) 

where A is a complex parameter. Autocovariance function becomes an analytically continued func- 
tional of the complex variable A: 

R~(u, A) — 2 S(f, A) cos(27r fu)df, (4) 
Jo 
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which coincides (due to the properties of analytically continued complex functions) exactly with 
the integral in ([!]) exept at the point A = 0. The functional (H) with the spectrum defined by eq. 
([|) is a table integral and can be easily calculated analitically. After calculation of the integral it 
is expanded in the Laurent series near the point A = 0. It yields 

R-(u,A) = -^Residue A=0 {iT(u, A) }+R~(u, A = 0) + ... . (5) 

The first term in the expansion is a simple pole with respect to A. The second term in the expansion 
is finite and gives exactly the autocovariance function given in eq. Regularization of the integral 
(H) means that we take only its finite part and abandon the singular term. Analytic continuation 
technique is powerful theoretical tool (Kopeikin 1997a) but it hardly can be used in numerical 
computations for which infrared cut-off technique is much better. 



2.2 Infrared Cut-off Technique 

The power spectrum S(f) of the noise is defined using the truncated Fourier transform with the 
lower cut-off frequency f — s. Namely, we require that the truncated cosine Fourier transform of 
S(f) must give the stationary part of the original autocovariance function (0) without any additional 
contributions. Let us postulate that the spectrum S(f) may be represented by the formula 



S(f) 



1 oo 

+ I>*(^ (2i) (/ 



k=0 



if/>e 
otherwise 



(6) 



where the spectral index of noise n = 1,2, ...,6, constant parameter h n is the strength of noise, 
quantities B^{e) are constant numerical coefficients being defined later, and 5^ k \f — e) denotes the 
n — th derivative with respect to / of the Dirac delta-function 5(f — e). The Dirac delta function 
is defined according to the condition (Korn & Korn 1968) 



f(x)S(x-X)dx 



0. 



if(x + o), 
If(x-o), 



if X < a, or X > b, 
if X = a, 
if X = 6, 



(7) 



t l[f(X-0) + f(X + 0)}, ifa<A<6, 



where f(x) is arbitrary function being such that unilateral limits f(X — 0) and f(X + 0) exist. 
Coefficients Bk(e) are determined by the condition that the cosine Fourier transform of S(f) gives 
stationary part of autocovariance function of corresponding low-frequency noise the model of which 
may be found in papers (Kopeikin 1997b, 1999). As an example of derivation of Bk{s) we determine 
the several first coefficients Bk(e) in the event of flicker noise in pulsar's rotational phase which 
has the spectral index n — 1. Coefficients B^{e) for noises having other spectral indices will be 
displayed in this section without proof which is rather straightforward. 

Stationary part of autocovariance function R~(t) of flicker noise in pulsar's phase is equal to 
/ii7r -1 log \u\. According to definition of the spectrum we should have 



POO 

R-(u) = 2 J S(f) cos(2tt fu)df. 
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Substituting S(f) from Eq. @ with n = 1 to right hand side of Eq. (Q) and taking integrals we get 



h 00 
R-(u) = -Ci(27re|u|) + h x cos(27rew) V(-l)^B 2fc (£)(2W) 2fc , (9) 

where Ci(x) is the cosine integral, and we have used the formula (Korn & Korn 1968) 

2 / s^if - e) cos(2tt fu)df = -7^- cos(2neu) = (-l) k (2nu) 2k cos(2neu) (10) 
J e ds 

Taylor expansion of the cosine integral and cos(2-7T£ti) in the right hand side of the expression @ 
with respect to small parameter e yields 



R-(u) = — [- log \u\-j- log(27T£)] + hA Bo(e) + (2neuf 



TT 



L _ l B ,(s) - B 2 (e) 
271 Z 



+ 0(e 4 ).(ll) 



Since R (u) must be equal to — hiir 1 log \u\, we find from Eq. (|TT| ) 



= 2 + clogM = f - f - ^ (12) 

7T 7T 477 Z7T Z7T 

For practical purposes it is enough to account for the coefficient B (e) only, since all other resid- 
ual terms appear being multiplied by the factor 2tt6u which is negligibly small under the usual 
circumstances because of smallness of the product eu. Thus, residual terms are not important as 
long as the accuracy of observations is not high enough. It is worth emphasizing that the residual 
terms under discussion are model dependent. Had we chosen another model for the spectrum having 
slightly another behavior as frequency approaches to zero the residual terms would look differently. 
These arguments permit to estimate how long we can observe a certain pulsar and process the data 
with one or another model of red noise spectrum. 

Proceeding in the same way for the set of other spectral indeces we obtain the following expres- 
sions for the power spectra of low-frequency noises: 

(1) Flicker noise in phase: 

S(f) = ^{j + 2 [7 + log(2vr5)] 5(f -e)} + 0(e 2 ). (13) 

(2) Random walk in phase: 

(3) Flicker noise in frequency: 

S(f) = ^- 3 {j 3 - ^S(f -e) + [log(27r £ ) + 7 -l] 6^(f - e)} + 0(s 2 ). (15) 

(4) Random walk in frequency: 
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W = ^-^C/-e)-^C/-e)} + 0( e) . (16) 



167T 4 I / 

(5) Flicker noise in frequency derivative: 



32tt 5 1/ 

(6) Random walk in frequency derivative: 

s <« = " ~ £) " ih^ f - £ >} + °< £ »- < 18 > 

Let us note that the coefficient B^{e) = in the expression fll8|) for spectrum of random walk in 
frequency derivative. 

The expressions given above indicate that there is a strong concentration of infinite energy of 
the noise at the lower cut-off frequency. As we have stressed already it is a specific feature of the 
chosen model of the spectrum which appears because we do not know a real behavior of spectrum 
while frequency is approaching to zero. Another remark is that a direct integration of any of the 
foregoing spectra with respect to frequency from / = e to infinity (which may be erroneosly treated 
as the energy being stored in TOA residuals) gives zero value which may look surprising. However, 
it is worth noting that the entire energy presents in TOA residuals can be calculated only after 
multiplication of the spectrum by the filter function (see section 6 below for more detail). Therefore, 
calculation of the total energy of residuals is more complicated and always gives a positive numerical 
value as it should be. Similar arguments can be used in calculating variances of fitting parameters. 
For example, calculation of variances of the first several spin-down parameters in frequency domain 
may give a negative numerical value of the variance (Kopeikin 1999) which is physically meaningful. 
The paradox is solved if we remember about contribution of the non-stationary part of noise which 
always makes variances of the parameters numerically positive (for more detail see (Kopeikin 1999)). 

Pulsar timing observations can be used for estimation of strength and spectrum of low frequency 
noise presents in TOA residuals. For this reason, development of practically useful estimators of 
spectrum of noise are required. We are not going to consider in the present paper the question about 
how to construct the best possible estimators. This subject has been enlightened by a number of 
other researches (see, for instance, the papers of Deeter & Boynton (1982), Deeter(1984), Taylor 
1991, Matsakis al. 1997). Our purpose is to study spectral dependence of TOA residuals and 
variances of fitting parameters which are used in real practice. In order to make clear what we are 
doing let us describe, first of all, the timing model we are dealing with. 



3 TIMING MODEL 

We consider a simplified, but still realistic model of arrival time measurements of pulses from a 
pulsar in a binary system. It is assumed that the orbit is circular, and the pulsar rotates around 
its own axis with angular frequency u p which slows down due to the electromagnetic (or whatever) 
energy losses. It is also taken into account that the orbital frequency of the binary system, n&, and 
its projected semimajor axis, x, have a secular drift caused by radial acceleration of the binary 
(Damour & Taylor 1991, Bell & Bailes 1996), its proper motion of in the sky (Kopeikin 1996), and 
emission of gravitational waves from the binary (Peters & Mathews 1963, Peters 1964) bringing 
about the gravitational radiation reaction force (Damour 1983a, Grishchuk & Kopeikin 1983). 

© 1999 RAS, MNRAS 000, 000-000 



Millisecond and Binary Pulsars as Nature 's Frequency Standards 7 



The moment T of emission of the jV-th pulsar's pulse relates to the moment t of its arrival mea- 
sured at the infinite electromagnetic frequency by the equations (Damour & Taylor 1992, Kopeikin 
1994, 1999): 

D[T + xsm(n b T + a)} = t + <p (t) + <ft(t), (19) 
t = r* + A c + A R0 + A wQ + A EQ + A 50 . (20) 
We use the following notations: 

• T - pulsar time scale, 

• t - barycentric time at the barycenter of the Solar system, 

• t* - topocentric time of observer, 

• Ac, A/jq, A nQ , Aeq, As q - clock and astrometric corrections (Taylor & Weisberg 1989, 
Doroshenko & Kopeikin 1990, 1995) which one assumes to be known precisely, 

• D - Doppler factor gradually changing due to the acceleration and proper motion of the binary 
system in the sky Q 

ma- initial (constant) orbital phase, 

• n b - orbital frequency {n b = 2ir/P h ), 

• i - angle of inclination of the orbit to the line of sight, 

• x - projected semimajor axis a p of the pulsar's orbit (x = a p smi/c), 

• c - speed of light, 

• (fo(t) - the gaussian noise of TOA measuring errors, 

• tpi(t) - low-frequency gaussian noise caused by the long-term instabilities of terrestrial clocks, 
effects in propagation of radio signals in the interstellar medium and stochastic background of 
primordial gravitational waves, etc. 

The rotational phase of the pulsar is given by the polynomial in time 

Af(t) = v p T+U p T 2 + U p T 3 + ± u p T 4 + v 9 T 5 + W (T) + O (T 6 ) , (21) 

where u p , u p , u p , etc. are pulsar's rotational frequency and its time derivatives all referred to the 
epoch T = 0, the term O (T 6 ) denotes high order derivatives of the rotational phase, and (f2{T) 
is the intrinsic pulsar timing noise in either rotational phase, frequency, or frequency derivative. 
Solving iteratively equation (|I9"D with respect to T and substituting T for the right hand side of 
equation (ETf) gives a relationship between two observable quantities M and t: 



Af{t) = Wo + vt + \v t 2 + I v t 3 + — v t 4 + 

2 o 24 

1 11 11 

v t 5 - v(x+ x t + - x t 2 + - x t 3 ) sin(cr + n b t + -h b t 2 + -h b t 3 ) + ue{t), (22) 



120 v 2 6 ' v 2 6 



e(t) = f Mt) + <Pi(t)+Mt), (23) 

where A/"o is the initial rotational phase of the pulsar (Ao — —vt Q ); u, is, is,... are the pulsar's 
rotational frequency and its time derivatives at the initial epoch i ; 

the projected 

semimajor axis of the orbit and its time derivatives at the epoch to; a, n b , n b , n b , ... are the pulsar's 
orbital initial phase, orbital frequency and its time derivatives at the epoch i - Timing model (|22|) 



D = — c a, where Vr and V are correspondingly the relative radial and total velocities of the binary system barycentre with 



respect to the barycentre of the Solar system 
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Table 1. List of the basic functions and parameters used in the fitting procedure. Spin parameters 
8Mo,8v, 8 u,8 v,8 u,8 v , fit rotational motion of the pulsar around its own axis. Keplerian parameters 
8x, 8a, 8ni, fit the Keplerian orbital motion of the pulsar about barycentre of the binary system. Post- 
Kcplcrian parameters 8 x, S x, 8 x ,8 n b ,8 n b fit small observable deviations of the pulsar's orbit from the 
Keplerian motion caused by the effects of General Relativity, radial acceleartion, and proper motion of 
barycentre of the binary system with respect to the observer 



Parameter Fitting Function 

fe = ^4r *2(t) = tt 
2rl b 

/34 = - i T 4r i/>4(t) = u 3 

^ = i±jNr <fe(i) = « 4 

fa 

(3j — — Sx sin a — Sax cos a ipj{t) — cos u 

/3g — — Sx cos a + Sax sin a i/>g (t) = sin u 

(3g — ^ — S x cos a A- Sn b x sin a^ ijjg (£) — u sin u 

010 — (^—Sx sin a — Sn^x cos a^j V'lO (*) = u cos u i 

/3l2 — — -j 1 ( —S x cos a - -f- S n fa a; sin a J ^12 (*) — ^ sinu 

/3i3 — — ^ir ( — S x cos cr + S rift x sin tr) ^13 (*) ~ 14 ^ sin u 

(3 14 — - 1 3 ^ — <5 a; sin a — S n b a; cos cr^ ^14 (*) — ^ 3 cos u 



accounts only for linear terms which is enough for implication of least square method of fitting 
parameters to the data. All non-linear residual terms of order x 2 , xe, e 2 , i>x, i)x, etc. are negligible 
and, for this reason, have been omitted from (p2|). 

We assume that all observations of the binary pulsar are of a similar quality and weight. Then 
one defines the timing residuals r(t) as a difference between the observed number of the pulse, 
J\f obs , and the number J\f(t,9), predicted on the ground of our best guess to the prior unknown 
parameters of timing model (|22|), divided by the pulsar's rotational frequency u, that is 



rM) = *»•-*(«■») , (24) 

where 9 = {9 a ,a = l,2,...k} denotes a set of k measured parameters [k = 14in the model fl22|)) 
which are shown in Table [l|. It is worth noting that hereafter we use for the reason of convinuence 
the time argument u = ribt, that is the current orbital phase. If a numerical value of the parameter 
9 a coincides with its true physical value 9 a , then the set of residuals would represent a physically 
meaningful noise e(t), i.e. 

r(t,9) = e(t). (25) 
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In practice, however, the true values of parameters are not attainable and we deal actually with 
their least square estimates 9*. Therefore, observed residuals are fitted to the expression which is a 
linear function of corrections to the estimates 9* a of a priori unknown true values of parameters 6 a . 
From a Taylor expansion of the timing model in equation (0), and the fact that r(t, 9) = e(t) one 
obtains 

r(t, 9*) = e(t) - ^.(i, 9*) + 0((3 2 a ), (26) 

a=l 

where the quantities (3 a = 89 a = 9* — 9 a are the corrections to the presently unknown true values 
of parameters, and ip a (t, 9*) = are basic fitting functions of the timing model. 

In the following it is more convenient to regard the increments (3 a as new parameters whose 
values are to be determined from the fitting procedure. The parameters /3 a and fitting functions are 
summarized in Table [I] with asterisks omitted and time t is replaced for convenience by the function 
u = n b t which is the current value of orbital phase. We restrict the model to 14 parameters since 
in practice only the first several parameters of the model are significant in fitting to the rotational 
and orbital phases over the available time span of observations. 

Let us introduce auxilary functions E a (t) defined according to the formula: 

E a (t) = EL-^t), (27) 

6=1 

where the matrix of information 

mN 

L ab (T) = J2MU)MU), (28) 

i=l 

the matrix L~ h l is its inverse, and T = NP^ is a total span of observational time. Functions S(t) are 
called (Deeter 1984) the dual ones to ip a {t) because of the cross-orthonormality condition is hold: 

mN 

Y^a(U)MU) = 5 ab . (29) 

i=l 

Now suppose that we measure m equally spaced and comparably accurate arrival times each orbit 
for a total of N orbital revolutions, so we have mN residuals r« = r(t»), i = 1, ...,mN. Standard 
least squares procedure (Bard 1974) gives the best fitting solution for estimates of the parameters 

A* 

mN 

f3 a (T) = Y^a(U)e(U), a = l, ...,14. (30) 
i=i 

Let the angular brackets denote an ensemble average over many realizations of the observational 
procedure. Hereafter, we assume that the ensemble average of the noise e(t) is equal to zero. Hence, 
the mean value of any parameter /3 a is equal to zero as well, i.e. 

< e (t) >= — ► < p a >= 0. (31) 
The covariance matrix M a & = < f3 a Pb > of the parameter estimates is now given by the expression 
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mNmN 

M ab (T) = Y / T, E a(t t )E b (t J )R(t l ,t j ), (32) 
i=ij=i 

where R(ti,tj) = < e(ti)e(tj) > is the autocovariance function of the stochastic process e(t). The 
covariance matrix is symmetric (that is, M a b = Mf, a ), elements of its main diagonal give variations 
(or dispersions) of measured parameters crg = M aa = < /3% >, and the off-diagonal terms represent 
the degree of statistic covariance (or correlation) between them. Covariance matrix consists of two 
additive components M^ b and M~ b describing correspondingly contributions from non-stationary 
and stationary parts of autocovariance function R(t{,tj). Explicit expressions for the matrix M a b 
can be found in the paper (Kopeikin 1988) wherein we have done all calculations in the time domain. 
Only M~ b admits transformation to the frequency domain which will be discussed in subsequent 
sections. 

Subtraction of the adopted model from the observational data leads to the residuals which 
are dominated by the random fluctuations only. An expression for the mean-square residuals after 
subtracting the best-fitting solution for the estimates (^) is given by the formula 

< r 2 (T) >= -VV^ i) i J )fl(i 1) { J ), (33) 



-r mNmN 

W >=—j vEE#Mi)*Mi). 



where the function 

14 

K(U, tj) = S tj - J2 E a{ti)Mtj), (34) 



a=l 



is called the filter function (Blandford et al. 1984). We have proved (Kopeikin 1999) that the post-fit 
residuals depend only on the stationary part of the noise 



-I mNmN i 14 14 

< At) >= = --.EE^ 1 

" tiv i=lj=l o=16=l 



mNmN 

i=ij=\ 



(35) 



For this reason, methods of spectral analysis in frequency domain can be applied for analyzing 
residuals without any restriction. Let us note that the explicit dependence of TOA residuals on the 
total span of observations contains in (Kopeikin 1988). 



4 FOURIER TRANSFORM OF FITTING FUNCTIONS 

We define the Fourier transform of the fitting functions ip a (t) as 

mN 

* a (/, m, N) = Y^Mtj) exp(-27ri/t i ), (36) 
i=i 

where / is the Fourier frequency measured in units being inversly proportional to units of mea- 
surement of time t. We measure time in units of orbital phase u = ntf, that is in radians. Then 
the frequency u = 2irf is dimensionless and measured in units of orbital frequency n&. One notes 
the Fourier transfrom of the fitting functions depends on the Fourier frequency /, total amount of 
orbital revolutions N, and frequency of observations m. 

When the total amount of observational points, mN, is large we can approximate the sum (J36f) 
by the integral (Kopeikin 1999) 
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(37) 



nN 



ttN 



ip a ( u ) exp(—itvu)du. 



(3* 



We note that ijj a (—ou) = ip*(u), where the asterisk denotes a complex conjugation. Replacing the 
sum over observational points by the integral with respect to time (or orbital phase) is equivalent 
to the case of continuous observations. 

The following formulae are also of use in practical computations: 



rnN 

2 / ip a (u) cos(uu)du, if index a = 1, 3, 5, 
Jo 



rirN 

—2i \ ipa(u) sm(uju)du, if index a = 2, 4, 6, 
*• Jo 



(39) 



These expressions shows that the fitting functions with odd indices are real and those with even 
ones do complex. 

Let us introduce notations - T = nN, z = uT. The Fourier transform of fitting functions takes 
the form: 



^13 (u 



2iT 2 Mz), 

2tT A Mz), 
2T 5 Mz), 

2tT e M*), 



= T 



r h (z + T)+Mz-T) 



tT[Mz + T)-Mz-T) , 
T 2 [fo( z -T)-fo(z + T) 
iT 2 [fa(z-T) + Mz + T) 
fo(z + T) + 4> 3 (z - T) 
'fo(z + T)-Mz-T) 
'Mz - T) - Mz + T) 



rj"i3 

iT 3 

rjn4 



iT 



l A (z - T) + 4> 4 (z + T) 



(40) 
(41) 
(42) 
(43) 
(44) 
(45) 
(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 

(53) 
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Functions 4> a {z) (a = 1,2,..., 6) have the following form: 

Hz) = , (54) 



>2\ 



«) = — -5n, (55) 



sin ^ 








z l 








2 cos 2 


2 sinz 






z 2 


z 3 ' 






3 sinz 


6 cosz 


6 sinz 




z 2 


z 3 






4 cos 2; 


12 sinz 


24 cos z 


24 sin z 


z 2 


z 3 


- 2 * + 


z° 


5 sinz 

z 2 


20 cos z 

z 3 


60 sin z 

+ + 


120 cos z 

z 5 



i 3 (z) = + — (56) 



*M = - — p — p-+— < 57 > 

7 . . sinz 4 cos z 12 sin z 24 cos z 24 sin z 

Hz) = + ^ = ^ + =— , (58) 



*M = —-^5 — + — + — i — • < 59 > 

Actually, it is more convenient to use in what follows the spherical Bessel functions j a (z) defined 
as (Korn & Korn 1968, section 21.8-8) 

^ = *'{-li)' S ¥' (-M.2,...) (60) 



Plots of the functions jo(z), ji(z),..., j$(z) are displayed in Figure (|Dl|). The functions have a dif- 
ferent behavior near the point z = 01, and then oscillate with monotonically decreasing amplitude. 
Asymptotic expansion of the sperical Bessel functions for large values of the variable z are given 
by the formula 

7TGf 



sin l z 

3a{z) « — * (61) 

z 

It is worth emphasizing that the maximal value of any of these functions can not be large than 1. 

Fitting functions (f> a {z) being expressed in terms of the spherical Bessel functions assume the 
form 

Hz) = Jo(z), (62) 
2 (z) = -h{z), (63) 

Hz) = ±M') ~ |*(*). (64) 
Hz) = -Iji(z) + Ij 3 (z), (65) 



t The function j a (z) = z + ° + 0(z a + 2 ) for z « 1 (a = 0, 1, 2, ...). 
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M z ) = ^jo(z) ~ jh{z) + 3^4(2), (66) 
M*) = -^h(z) + ^h(z)--^j 5 (z). (67) 

5 FOURIER TRANSFORM OF THE COVARIANCE MATRIX 

In order to calculate the covariance matrix we need to know the Fourier transform of the dual 
functions E a (t). The transform is defined in accordance with definition (|27|) of the dual functions 
and takes the form 

3 a (/, m, N) = EL-^ C (/, m, N), (68) 

c=l 

and the cross-orthonormal condition in frequency domain is given by the integral 

S (/, m, N)* b (f, m, N)df = -5 ab . (69) 
2 

In the limit of continuous observations it is convinuent to introduce matrix C ab = —L ab instead 
of the matrix L ab . Explicit expression for the matrix C ab is given by the integral: 

/wN 
ipa(u)ipb(u)du, (70) 
-TrAf 

and the result of calculation of the integral is given in the paper (Kopeikin 1999, Tables 5,6). Then 
we have L~ b = ^C~ b and the dual function H(/) can be recast as 

^a{f,N) = Y.C- b l Mf,N). (71) 
b=i 



Hence, comparing the eq. (71) with (|68| ) one concludes that in the limit of continuous observations 
the Fourier transfrom of the dual functions depend only on the Fourier frequency and total amount 
of orbital revolutions as it was expected. It is more insightful to express the dual functions ([7T| ) in 
terms of the spherical Bessel functions (§^). The expressions obtained are rather unwieldy and, for 
this reason they are given in Appendix A. Plots of the Fourier transform of the spherical Bessel 
functions are given in Appendix D. 

Making use of definition of the Fourier transforms of stationary part of autocovariance function 



(H) and the dual functions ( 68|) we obtain the Fourier transform of stationary part of the covariance 
matrix M~ b (m, N) 



M~ b = ^ S(f)H ab (f,™,N)df, (72) 
where H ab (f, m, N) is the transfer function given by the expression 

H ab (f, m, N) = E a (f, m, N)E* b (f, m, N) + E* a (f, m, N)E b (f, m, N), (73) 

and asterisk denotes a complex conjugation. For numerical computations of M~ b the next formula 
can be used in practice 

M- b (m,N) = J^-^H ab {f, m) N)f- n df+ (74) 
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-fl r 

2 



B (e)H ab (e,m, N) + e 2 B 2 (e)H%> (e,m, N) + e 4 S 4 ( £ )^(£, m, iV) + ... 



ab 



where derivatives of H a b are taken with respect to the Fourier frequency, ellipses denote terms of 
negligible influence on the result of the computation, and A is the upper cut-off frequency arising 
from the sampling theorem and inversly proportional to the minimal time between subsequent 
observational sessions. 



6 FOURIER TRANSFORM OF RESIDUALS 

Fourier transform of timing residuals is obtained from eqs. 



and (|38l) . This yields: 



S(f)K(f,m,N)df, 



(75) 



where K (/, m, N) is the Fourier transform of the filter function (|34j ) 

K(f, m, N) — 1 — [S a (/, m, N)** a {f, m, N) + E* a (f, m, N)* a (J, m, N) 



(76) 



In the limit of continuous observations there is no dependence on the frequency of observations, m, 
so that one obtains 

i 14 

K(f, N) — 1— —Y p a (f, N)i,* a (f, N) + E* a (f, N)Uf, N)] . (77) 

Plot of the Fourier transform (|77| ) of the filter function K(f) is shown in Appendix E for different 
amount of orbital revolutions N. It is approximately equal to 1 until frequency is higher than 1/T, 
and rapidly decreases its amplitude as frequency approaches to zero. We also note that the curve 
of the Fourier transform clearly shows the additional dip near the orbital frequency. The dips near 
zero and orbital frequencies are getting narrower as a number of observational points increases. 



7 SPECTRAL SENSITIVITY OF TIMING OBSERVATIONS OF MILLISECOND AND BINARY 
PULSARS 

Analytical expressions and graphical representations of Fourier transforms of dual functions and 
timing residuals help us to understand in more detail spectral sensitivity of single and binary pulsars 
to different frequency bands in spectral decomposition of noise. First of all, let consider behavior 
of Fourier transform of fitting functions near zero and orbital frequencies. 

It is easy to confirm after making use of Taylor expansion of exponential function in ( |38l) near 
uj ~ that 



C al - WCat + ±(jJ*C a5 + U 6 Pa 



if a=l,3,5,... 



-ujC a2 + |^ 3 C a4 - j^uj 5 C a6 + uj 7 Po) , if a=2,4,6,..., 



(75 



where p a is a residual term depending only on the total amount of orbital revolutions, N. Taylor 
expansion of Fourier transform of fitting functions near the orbital frequency yields: 



C a7 + (u- 1)C„9 
-C a8 + (u- l)C o . 10 



C a .l 



2 ^ a - a 6 

<^a.l2 



C a .!3 + l) 4 <?a, 



C a .u + (w - l) 4 q a 



if a=l,3,.. 



if a=2,4, 



(79) 
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where q a is a residual term depending only on the total amount of orbital revolutions, N (let 
us remind that frequency is measured in units of orbital frequency n b ). Applying to Eqs. (|78|)- 
([75|) definition of the dual functions ( [71]) in the limit of continuos observations yields asymptotic 
behavior of the dual functions 

r 5 al - \ui 2 5 a z + ±Lj A 5 a5 + u 6 P a , if a=l,3,5,... 

E a (u) = l (80) 
[ i [-ujd a2 + |w 3 5 a4 - j^uJ 5 d~ a6 + uj 7 P a ) , if a=2,4,6,..., 

near zero frequency, and 

5 a7 + {u- l)5 a9 - ^C a . u - ^5 a . 13 + l) 4 Q a , if a=l,3,... 



-6 a8 + l)5 a . 10 + ^p-5 aS 2 - (j ^ L S a . u + (w - l) 4 Qj , if a=2,4,..., 



2 

near the orbital frequency, where P a and Q a are residual terms. Table 2 shows the asymtotic 
behavior of the residual terms of the dual functions. 

Now we can study asymptotic behaviour of filter function K(f) defined by eq. (|77|) . Taking into 
account the fact that J2l=i ^ab^bc = o~ac - the unit matrix, we obtain 

B • u 6 , when uj — > 

h-(f) = { (82) 
Bi • (a? — l) 4 , when u — > 1 

where So and -Bi are numerical constants. Such specific behavior of the filter function K(f) sig- 
nificantly reduces the amount of detected red noise below the cut-off frequency f c ~ « C T _1 and in 
the frequency band 1 — a^T -1 < / < 1 + a^T -1 lying near the orbital frequency. Here constant 
coefficents a c and a b can be determined by means of comparision of calculations of mean value of 
timing residuals in time and frequency domains (Kopeikin 1997a). Low frequencies of noise power 
spectrum are fitted away by the polynomial fit for the spin-down parameters of the observed pulsar. 
Frequencies being close to the orbital one are fitted away by the fit for orbital parameters of the 
pulsar. Amount of noise power remained in timing residuals after completion of fitting procedure 
is estimated by the expression 

<r*>=2/ ¥ W + ^W^^+O^). (83) 

The second term in the right hand side of eq. (|83"D shows amount of noise absorbed by fitting orbital 
parameters. It is negligibly small comparatively with first term in the right hand side and can be 
not taken into account in practice. Hence, we declare that the post-fit timing residuals can be used 
for estimation of amount of red noise and its spectrum in frequency band just from a c T _1 up to 
infinity, irrespectively of whether the pulsar is binary or not. This reasoning puts on firm ground the 
estimates of spectral window of timing observations and cosmological parameter Q g , characterizing 
energy density of stochastic gravitational waves in early universe, made by Kaspi et al. (1994) and 
Camilo et al. (1994) using observations of binary pulsars PSR B1855+09 and PSR J1713+0747 
respectively. 

Analysis of spectral sensitivity of estimates of variances of spin-down and orbital parameters 
is more cumbersome. We are interested in which frequencies give the biggest contribution to the 
variances. This is important to know, for example, in the event of using variances of certain orbital 
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Table 2. Asymptotic behavior of residual terms of the dual functions H a near zero and orbital 
frequencies. Constant h = cosT = (— 1) . 



Dual function Residual term P a Residual term Q a 



T 6 12T 2 h 
E 56 



rT 6 -iT 2 h 



_J_T 4 -iZh 

1584 4 



. 1 T 4 --h 

6864 4 

_^j_ T 2 _45 h 

528 1 8 



1 T 2 33 h 

3120 40 

1 ^4}, 1 ^p4 

165 280 



-T 6 h -J_T 4 

5 280 



,^_T 4 h 3_T 2 

693 28 



-^_T 4 h J-T 2 

273 28 



19 rp2u 1 rp2 

693 28 



-J— T 4 h 3_T 2 

9009 28 

-J_T 2 h i 

297 12 



31 T 2 h _J_ 

3861 12 



parameters for setting the fundamental upper limit on Q g (Kopeikin 1997a, Kopeikin & Wex 1999). 
Analitic calculations reveal the leading terms in asymptotic expansions of the dual functions near 
zero frequencies: 

Squares of the dual functions £ a appearing in eqs. ( |B1| )-( [BT4|) and eqs. (|Cl|)- ([CT^ ) is rather 
complicated. Their behavior is periodic with bumps both near zero and orbital frequencies with 
rapidly decaying occilating wings far outside these frequencies. We evaluated that the frequency 
bump for the spin-down parameters near zero frequency is bigger than that near the orbital one. On 
the other hand, the situation for orbital parameters is just opposite. Analitic behavior of the dual 
functions shows that there are two spectral windows in which they are the most sensitive to the 
stochastic noise. One is located near zero frequency and the second one lies near the orbital one. We 
can bound these windows by two frequency intervals (0,^) and (1 — 1 + respectively where 
constant coefficients a, ct_, a + can be calculated by means of comparison of results of calculation 
of variances in time and frequency domains. It allows easily to calculate contributions from the 
foregoing frequency bands to numerical values of variances of fitting parameters and determine 
which of these frequency bands makes bigger deposit. Variances of fitting parameters depend only 
on the total span, T, of observations (Kopeikin 1997b). Comparing dependence of the variances 
on T calculated in time domain with that calculated in two frequency windows reveal the relative 
importance of different frequencies. In order to do this we have to compare magnitude of two 
integrals 
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Table 3. Comparative contribution of different frequency bands to variances of spin-down (a = 
1,2,..., 6) and orbital (a = 7, 8, ...,14) parameters f} a - Number n = 1,2,..., 6 denotes the spectral 
index of corresponding red noise. Time dependence of all variances completely coincides with that 
which was obtained by calculations in time domain as given in (Kopcikin 1999). 



Variance of parameter Contribution of integral Ii Contribution of integral I2 



A ^T"" 1 ~t- 3 

PI 

T 2 ~ T n " 3 ~ T" 5 

P2 

t% ~ T n_5 ~ T -7 

P3 



P5 

<j\ ~ rpn-11 ~ rp-13 

P6 

0-2 ^T™" 5 ~ T" 1 

P7 

PS 

ci ~ T' 1-5 ~ T~ 3 

P9 

t 2 , ~ T"" 7 ~ T" 3 

P10 

T 2 ^ T" -9 ~ T" 5 

P11 

T 2 ~ T n " 7 ~ T" 5 

P12 

T 2 ^ T "-9 ^ T -7 

P13 

T 2 ^ rpn-11 _ T -7 

P14 



and 



(27T/) 



fc=0 



34) 



I, 



(27T) 



\z a (f)\ 2 r n df, 



(85) 



The first integral describes contribution of low frequencies to parameter's variances. The second 
integral gives contribution to parameter's variances from the frequencies lying near the orbital 
one. It is worth emphasizing that terms with delta function and its derivatives in (|84| ) bring on 
mutual cancellation of all terms depending on cut-off frequency e and diverging as e goes to zero. 
Such cancellation has been expected since we modified the spectrum of red noise so that to avoid 
appearance of all divergent terms which have no physical meaning. We don't give here results of 
calculation of numerical values of parameters a, a_, a + because they are not so important for 
making conclusions. Asymptotic behavior of the dual functions near zero and orbital frequency is 
enough to see which frequncy band is the most important for giving contribution to corresponding 
integrals and parameter's variances. Time dependence of two integrals is shown in Table ([3]) up 
to not so important constant. Behavior of variances of the first three spin-down parameters is 
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not interesting for they are contaminated by the presence of the non-stationary part of red noise 
(Kopeikin 1997b). For the rest of spin-down parameters one observes that contribution of noise 
energy from low frequencies to variances of the parameters is dominating. However, the situation is 
not so simple in the event of orbital parameters. One can see that in the event of red noise having 
spectral index n < 4 contribution of the noise energy from the orbital frequency span (1 — ^,1 + ^) 
can be equal to or even bigger than that from the low frequency band. Only when the spectral index 
of noise n > 5 contribution of the noise energy of low frequencies to variances of orbital parameters 
begins to dominate. It is worth noting that the timing noise with spectral index n = 5 is produced 
by cosmological gravitational wave background. The fact that for this noise low-frequencies give the 
main contribution to variances of orbital parameters confirms our early statement (Kopeikin 1997a) 
that measurement of variances of orbital parameters showing secular evolution tests the ultra-low 
frequency band of cosmological gravitational wave background. Hence, these variances can be used 
for setting upper limit on the cosmological parameter Q g in this frequency range in contrast to 
timing residuals which test only low-frequency band of the background noise. (Kopeikin 1997a). 
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APPENDIX A: EXPLICIT EXPRESSIONS FOR THE DUAL FUNCTIONS 

In this appendix we give explicit expressions for the dual functions. Using definition ([FT]) and 
elements of inverse matrix from the paper (Kopeikin 1999, Tables 5 and 6) one obtains 



_ 5 27 

= j (z) + -j 2 (z) + -^h{z) + 

15h f 



-|3 [ 3l (z + T) - 3l (z — T)] — 7 [j 3 (z + T) - j 3 (z - T)] 

45h f 1 

•— 13 \j (z + T) +j (z - T)] - 20 [j 2 (z + T) + j 2 (z - T)]| , 



(Al) 
(A2) 
(A3) 



1^ 



^{z 



-E 4 (z) 
i 



3 
T 



3i{z) + -j 3 (z) + yJs(^) 



+ 



105h f 



{j Q (z + T) - 3o (z - T) - 5 [j 2 (z + T) - j 2 (z - T)] } + 



105h f 



— 151 \j x {z + T) + 3l (z - T)] - 154 [j 3 (z + T) + j 3 (z - T)] 



15 



2T 2 

105h f 



^{3 \ji(z + T) - Uz - T)] - 7 \j 3 (z + T) - j 3 (z - T)]| + 

1 D5h f 1 

^{7 [*>(* + T) +j (z - T)] - 50 [j 2 (z + T) + j 2 (z - T)]} , 



35 
2T3 



'• / N 11./ \ 



315h f 



j>(* + T) - j (* - T) - 5 \j 2 (z + T) - j 2 (* - T)] 
[Ji(z + T)+Uz - T)] - 28 [fcC* + T) + j 3 (z - T)] 



(A4) 
(A5) 
(A6) 

(A7) 
(A8) 
(A9) 

(A10) 
(All) 
(A12) 



© 1999 RAS, MNRAS 000, 000-000 



20 Sergei M. Kopeikin and Vladimir A. Potapov 
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SbC*) = W4 H{z) (A13) 

+ ^{3 \ji(z + T) - 3l (z — T)] — 7 + T) - j 3 (z - T)]} (A14) 

1 575h r 1 

~ 8Ter{[^ + T) + jo(* - T)] - 8 + T) + j 2 (z - T)]} , (A15) 

tS 6 (z) = Jpj 5 (*) + (A16) 

8T«r{ + T) - j (z - T) - 5 [j 2 (, + T) - - T)] } + (A17) 

^?{l3 [Ji(z + T) + - T)] - 42 b' 3 (^ + T) + h(z - T)]} , (A18) 

S 7 (z) = Jo ( 2 + T)+jo(^-T) + ^[j 2 (^ + T)+j 2 (^-T)] (A19) 

-^{3 LM* + T) - Mz — T)] — 7 [j 3 (z + T) - - T)]| (A20) 

[*»(*) - 6j 4 W] , (A21) 

-S 8 (z) = j ( z + T)-j (z-T) + ^\j 2 (z + T)-j 2 (z-T)}+ (A22) 

^{3 [ 3l (z + T) + - T)] - 7 [j 3 (z + T) + - T)]} + (A23) 
3h 

Y [3 jl (z)-7Uz) + Uj 5 (z)} , (A24) 

S»(*) = Ilii^ + ^-A^-^ + ^s^ + T)-^-^]} (A25) 

-^{bo(^ + T) + j>(* - T)] - 5 [j 2 (z + T) + j 2 (z - T)]} (A26) 
1 5h 

-^[>o(*)-5j 2 (z) + 9j 4 (*)] , (A27) 

tS 10 (*) = —{bi^ + ^+^-T^ + ^^ + Tj+^-T)]} (A28) 

-^{j'o(^ + T) - - T) - 5 [j 2 (z + T) - j 2 (z - T)] + (A29) 
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1 5h 

— [-m(z) + 77 33 (z) - 220j 5 (z)] , (A30) 

SnW = -^- 2 \j2(z + T)+j 2 (z-T)} (A31) 

+^{ 3 + T ) " " T )] " 7 + T) - j 3 (z - T)]} + (A32) 
1 5h 

— [2 Jo (,)+5j 2 (,)-72 j4 (^)] , (A33) 

t3 12 (s) = -^- 2 [h(z + T)-j 2 (z-T)} (A34) 
-^{3 + T) +MZ - T)] - 7 + T) +.73(2 - T)]} + (A35) 

1 ^ih 

[3jiW - 7i 3 (^) + Uj 5 (z)] , (A36) 

35 

= -^L73(^ + T)-^-T)]+ (A37) 

j^{jo(* + T) + j (* - T) - 5 [j^ + T) + j 2 (z — T)]| + (A38) 
35h 

— [jo ( z )-5j 2 (z) + 9j 4 (z)} , (A39) 
1 35 

tSuW = ^l>s(* + T)+j s (z-T)] + (A40) 

J^oC* + T) - j (* - T) - 5 y 2 (z + T) - - T)]} + (A41) 
1 05h 

— [4^-21^) +66j 5 (z)] • (A42) 

APPENDIX B: ASYMPTOTIC BEHAVIOR OF THE DUAL FUNCTIONS NEAR ZERO FREQUENCY 

In this appendix we give asymptotic behavior of the dual functions near zero frequency. They are 
as follows: 

Si(/) = L(z), £i(z)=j (z) + lj 2 (z) + ™u(z), (Bl) 
1 ~ 3 ~ ~ 7 55 

-S 2 (/) = Mz) = Mz) + -j 3 (z) + -j 5 (z) , (B2) 

S 3 (/) = Uz)=J2(z) + 9 -Uz), (B3) 
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\Uf) = §- 3 Uz), Uz)=h(z) + ^-j 5 (z), (B4) 
SI ^ ~ 

5 5 (/) = g^&W, &(*)=*(*), (B5) 
1 ~ fiQS ~ 

-2 6 (7) = -g^&W, (B6) 

S 7 (/) = -^£ 7 (*), e»=J 2 W-6j 4 W-^sinz, (B7) 

1~ 9h~ ~ 7 11 1 

tS 8 (/) = —6^), M^ii^-g^ + yjs^-gSinz, (B8) 
15h ~ 

S 9 (/) = -—£ 9 (z), ^(z)=j (z)-5j 2 (z) + 9u(z)-coaz, (B9) 

U 270h~ , . ~ 77 110 5. 1 

T=io(/) = — ^-CioW , ^ 10 {z)=j 1 (z)-— j 3 (z) + —j 5 (z)-—sm z-—z cos z , (BIO) 

30h ~ 5 1 

Sn(/) = -^4-^11(2), £n(z) =j (z) + -j 2 (z) -3Qj 4 (z) - -zsinz-cosz , (Bll) 

jS 12 (/) = -^£i 2 (*), i 12 (z)=j 1 (z)- 7 -j 3 (z) + ^j 5 (z)-^smz , (B12) 

Si 3 (/) = -^M*), £isW=6W, (B13) 

1~ 420h~ . . - . . . . 63 , . 33 , . 1 1 

t=i 4 (/) = -7^14(2), £14(3) = M*) - + -jM z ) ~ l smz ~ w zcosz ' ( B14 ) 



where h = (— 1) 



APPENDIX C: ASYMPTOTIC BEHAVIOR OF THE DUAL FUNCTIONS NEAR ORBITAL 
FREQUENCY 

Leading terms in asymptotic expansions of the dual functions S a (z) near orbital frequency are as 



follows: 

45h~ ~ 7 1 

= — , Zi(y) = ji(y) - ^k\y) - 3 siny , (ci) 

^S 2 (/) = -^| 2 (y), |2(y)=Jo(y)-5j 2 (y)-oosy, (C2) 

5 3 (/) = -^la(y), fs(v) = £i(v) (C3) 
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-Mf) = -^ld(y), Uv) = Uy), (C4) 

H 5 (/) = ^f 5 (y), f 5 (y)=fi(y), (C5) 

7 Se(/) = ^e(y), &(v) = 6(i/), (C6) 

Sr(/) = fr(y), | 7 (y)=io(y) + ^2(y), (C7) 

-S 8 (/) = f g (y), | 8 (y) = -| 7 (y) , (C8) 
% 

s 9 (/) = -|f 9 (y), lo(y)=ji(y) + ^s(y), (C9) 

tS w (/) = Uw(y), Uy) = 6(y), (Cio) 

v\ _ 

Sn(/) = -^n(y), 6i(y)=j 2 (y), (en) 

t3 u (/) = ^li 2 (y), li 2 (y) = -In(y) , (C12) 

Sis(/) = 2fs&3(v), ei 3 (y)=j 3 (y), (Ci3) 

^s 14 (/) = ^fu(y), Ii4(y) = li 3 (y) . (C14) 

where y = z — T. 

APPENDIX D: PLOTS OF THE FOURIER TRANSFORM OF SPHERICAL BESSEL FUNCTIONS 
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Figure Dl. Plots of the Fourier transforms of the sperical Bessel functions jo(z), ji(z),. ■ . , j§{z) in terms of the variable z = uiT. Ciclic 
frequency ui is measured in units of orbital frequency n\y. Amplitude of the transform is normalized to unity. 



APPENDIX E: PLOTS OF THE FOURIER TRANSFORM OF AUTOCOVARIANCE FUNCTION 
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Figure El. Plot of the Fourier transform of the filter function of timing residuals for the amount of orbital revolutions N = 4. Frequency 
is measured in units of orbital frequency n^. Amplitude of the transform has been normalized to unity 
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Figure E2. Plot of the Fourier transform of the filter function of timing residuals for the amount of orbital revolutions N = 10. Frequency 
is measured in units of orbital frequency n^. Amplitude of the transform has been normalized to unity. 
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Figure E3. Plot of the Fourier transform of the filter function of timing residuals for the amount of orbital revolutions N = 30. Frequency 
is measured in units of orbital frequency rtj. Amplitude of the transform has been normalized to unity. 
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